<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><meta name="generator" content="MATLAB 2021a"><title>Numerical properties of a reconstruction</title><style type="text/css">.rtcContent { padding: 30px; } .S0 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 28.8px; min-height: 0px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 24px; font-weight: normal; text-align: left;  }
.S1 { margin: 2px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S2 { margin: 20px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.S3 { margin: 3px 10px 5px 4px; padding: 0px; line-height: 20px; min-height: 0px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 20px; font-weight: bold; text-align: left;  }
.CodeBlock { background-color: #F7F7F7; margin: 10px 0 10px 0;}
.S4 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S5 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 4px 4px 0px 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S6 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 0px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S7 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px 0px 4px 4px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S8 { margin: 10px 10px 9px 4px; padding: 0px; line-height: 21px; min-height: 0px; white-space: pre-wrap; color: rgb(0, 0, 0); font-family: Helvetica, Arial, sans-serif; font-style: normal; font-size: 14px; font-weight: normal; text-align: left;  }
.S9 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 0px none rgb(0, 0, 0); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 0px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S10 { color: rgb(64, 64, 64); padding: 10px 0px 6px 17px; background: rgb(255, 255, 255) none repeat scroll 0% 0% / auto padding-box border-box; font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px; overflow-x: hidden; line-height: 17.234px;  }
.S11 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 0px none rgb(0, 0, 0); border-radius: 0px; padding: 6px 45px 0px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S12 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S13 { border-left: 1px solid rgb(233, 233, 233); border-right: 1px solid rgb(233, 233, 233); border-top: 1px solid rgb(233, 233, 233); border-bottom: 1px solid rgb(233, 233, 233); border-radius: 4px 4px 0px 0px; padding: 6px 45px 4px 13px; line-height: 17.234px; min-height: 18px; white-space: nowrap; color: rgb(0, 0, 0); font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 14px;  }
.S14 { margin: 10px 0px 20px; padding-left: 0px; font-family: Helvetica, Arial, sans-serif; font-size: 14px;  }
.S15 { margin-left: 56px; line-height: 21px; min-height: 0px; text-align: left; white-space: pre-wrap;  }
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsErrorElement {    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.embeddedOutputsErrorElement.rightPaneElement {}
/* Styling that is common to warnings and errors is in diagnosticOutput.css */.embeddedOutputsWarningElement{    min-height: 18px;    max-height: 250px;    overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsWarningElement.rightPaneElement {}
/* Copyright 2015-2019 The MathWorks, Inc. *//* In this file, styles are not scoped to rtcContainer since they could be in the Dojo Tooltip */.diagnosticMessage-wrapper {    font-family: Menlo, Monaco, Consolas, "Courier New", monospace;    font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {    color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {    color: rgb(255,100,0);    text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {    color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {    color: rgb(230,0,0);    text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart,.diagnosticMessage-wrapper .diagnosticMessage-causePart {    white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {    white-space: pre;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {    white-space: pre;    word-wrap:  initial;    min-height: 18px;    max-height: 250px;    overflow: auto;}
.textElement,.rtcDataTipElement .textElement {    padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement,.embeddedOutputsVariableStringElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement.rightPaneElement,.embeddedOutputsVariableStringElement.rightPaneElement {    min-height: 16px;}
.rightPaneElement .textElement {    padding-top: 2px;    padding-left: 9px;}</style></head><body><div class = rtcContent><h1  class = 'S0'><span style=' font-weight: bold;'>Numerical properties of a reconstruction</span></h1><div  class = 'S1'><span style=' font-weight: bold;'>Authors: Laurent Heirendt, Ronan M.T. Fleming, Luxembourg Centre for Systems Biomedicine</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Reviewers: Sylvain Arreckx, Thomas Pfau, and Catherine Fleming,  Luxembourg Centre for Systems Biomedicine</span></div><h2  class = 'S2'><span>Introduction</span></h2><div  class = 'S1'><span>During this tutorial, you will learn how to determine and explore the numerical properties of a stoichiometric matrix. The numerical properties are key to analyzing the metabolic reconstruction at hand, to select the appropriate solver, or to determine incoherencies in the network. </span></div><div  class = 'S1'><span>You will also be able to learn more about the definitions of the various numerical characteristics. This tutorial is particularly useful when you have a multi-scale model and are facing numerical issues when performing flux balance analysis or any other variants of FBA.</span></div><h2  class = 'S3'><span>EQUIPMENT SETUP</span></h2><h2  class = 'S2'><span style=' font-weight: bold;'>Initialize the COBRA Toolbox.</span></h2><div  class = 'S1'><span>Please ensure that The COBRA Toolbox has been properly installed, and initialized using the </span><span style=' font-family: monospace;'>initCobraToolbox</span><span> function.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span >initCobraToolbox(false) </span><span style="color: rgb(2, 128, 9);">% false, as we don't want to update</span></span></div></div></div><h2  class = 'S3'><span>PROCEDURE </span></h2><div  class = 'S1'><span>TIMING: 5 seconds - several hours (depending on the model size)</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Define the name of the model</span></div><div  class = 'S1'><span>Throughout this tutorial, we will use the </span><span style=' font-style: italic;'>E.coli core</span><span> model [2]. It is generally good practice to define the name of the file that contains the model, the name of the model structure, and the name of the stoichiometric matrix, as separate variables. Therefore, we propose that within the </span><span style=' font-family: monospace;'>modelFile</span><span>, there is a structure named </span><span style=' font-family: monospace;'>modelName</span><span> with a field </span><span style=' font-family: monospace;'>matrixName</span><span> that contains the stoichiometric matrix </span><span style=' font-family: monospace;'>S</span><span> (or </span><span style=' font-family: monospace;'>A</span><span>).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% If this is a distributed model, get the folder for the model. Commonly if you use your own model, this is unnecessary</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% As the file would be in the current folder. But for this tutorial we need to make sure, that the right model is used, </span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% and that no other model with the same name is above it in the path.</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >modelFolder = getDistributedModelFolder(</span><span style="color: rgb(170, 4, 249);">'ecoli_core_model.mat'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% define the filename of the model</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >modelFile = [modelFolder filesep </span><span style="color: rgb(170, 4, 249);">'ecoli_core_model.mat'</span><span >];</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% define the name of model structure</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >modelName = </span><span style="color: rgb(170, 4, 249);">'model'</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% define the fieldname of the stoichiometric matrix</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >matrixName = </span><span style="color: rgb(170, 4, 249);">'S'</span><span >;</span></span></div></div></div><div  class = 'S1'><span style=' font-weight: bold;'>Load the stoichiometric matrix</span></div><div  class = 'S1'><span>In order to use the model, we need to read the </span><span style=' font-family: monospace;'>modelFile</span><span> that contains a COBRA model structure  </span><span style=' font-family: monospace;'>modelName</span><span>:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% load the modelName structure from the modelFile</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >model = readCbModel(modelFile, </span><span style="color: rgb(170, 4, 249);">'modelName'</span><span >,</span><span style="color: rgb(170, 4, 249);">'model'</span><span >);</span></span></div></div></div><div  class = 'S8'><span>Some models contain stoichiometric matrices with a different name (commonly coupled models). By default, the stoichiometric matrix is denoted </span><span style=' font-family: monospace;'>S</span><span>.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% select the matrix</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >S = model.S;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">if </span><span >isfield(model, matrixName) == 1 &amp;&amp; strcmp(matrixName, </span><span style="color: rgb(170, 4, 249);">'A'</span><span >) == 1</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    S = model.A;</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div></div><div  class = 'S1'><span style=' font-weight: bold;'>Basic numerical characteristics</span></div><div  class = 'S1'><span>The </span><span style=' font-weight: bold;'>number of elements</span><span> represents the total number of entries in the stoichiometric matrix (including zero elements). This number is equivalent to the product of the number of reactions and the number of metabolites.</span></div><div  class = 'S1'><span>The number of rows represents the </span><span style=' font-weight: bold;'>number of metabolites</span><span> in the metabolic network. The number of columns corresponds to the </span><span style=' font-weight: bold;'>number of biochemical reactions</span><span> in the network.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the number of reactions and metabolites in S</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >[nMets, nRxns] = size(S)</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>nMets = 72</div><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>nRxns = 95</div></div></div><div class="inlineWrapper"><div  class = 'S11'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the number of elements in S</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >nElem = numel(S)  </span><span style="color: rgb(2, 128, 9);">% Nmets * Nrxns</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>nElem = 6840</div></div></div></div><div  class = 'S8'><span>The total number of nonzero elements corresponds to the total number of nonzero entries in the stoichiometric matrix (excluding zero elements).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the number of nonzero elements in S</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >nNz = nnz(S)</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>nNz = 360</div></div></div></div><div  class = 'S1'><span style=' font-weight: bold;'>Sparsity and Density</span></div><div  class = 'S1'><span>The </span><span style=' font-weight: bold;'>sparsity ratio</span><span> corresponds to a </span><span>ratio of the number of zero elements and the total number of elements. The sparser a stoichiometric matrix, the fewer metabolites participate in each reaction. The sparsity ratio is particularly useful to compare models by how many metabolites participate in each reaction.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the sparsity ratio of S (in percent)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >sparsityRatio = (1 - nNz / nElem) * 100.0  </span><span style="color: rgb(2, 128, 9);">% [%]</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>sparsityRatio = 94.7368</div></div></div></div><div  class = 'S8'><span>Similarly, the </span><span style=' font-weight: bold;'>complementary sparsity ratio</span><span> is calculated as the difference of 100 and the sparsity ratio expressed in percent, and therefore, is a ratio of the number of nonzero elements and the total number of elements.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the complementary sparsity ratio (in percent)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >compSparsityRatio = 100.0 - sparsityRatio  </span><span style="color: rgb(2, 128, 9);">% [%]</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>compSparsityRatio = 5.2632</div></div></div></div><div  class = 'S8'><span>The</span><span style=' font-weight: bold;'> average column density </span><span>corresponds to a ratio of the number of nonzero elements in each column (i.e. reaction) and the total number of metabolites. The average column density corresponds to the arithmetic average of all the column densities (sum of all the reaction densities divided by the number of reactions).</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% add the number of non-zeros in each column (reaction)</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >colDensityAv = 0;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">for </span><span >j = 1:nRxns</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    colDensityAv = colDensityAv + nnz(S(:, j));</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% calculate the arithmetic average number of non-zeros in each column</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >colDensityAv = colDensityAv / nRxns   </span><span style="color: rgb(2, 128, 9);">% [-]</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>colDensityAv = 3.7895</div></div></div></div><div  class = 'S8'><span>The average column density provides a measure of how many stoichiometric coefficients participate in each biochemical reaction on average.</span></div><div  class = 'S1'><span>The </span><span style=' font-weight: bold;'>relative column density</span><span> corresponds to the ratio of the number of nonzero elements in each column and the total number of metabolites. The relative column density corresponds to the average column density divided by the total number of metabolites (expressed in percent). The relative column density may also be expressed as parts-per-million [ppm] for large-scale or huge-scale models. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the density proportional to the length of the column</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >colDensityRel = colDensityAv / nMets * 100  </span><span style="color: rgb(2, 128, 9);">% [%]</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>colDensityRel = 5.2632</div></div></div></div><div  class = 'S8'><span>The relative column density indicates how many metabolites are being used on average in each reaction relative to the total number of metabolites in the metabolic network.</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Sparsity Pattern (spy plot)</span></div><div  class = 'S1'><span>The visualisation of the sparsity pattern is useful to explore the matrix, spot inconsistencies, or identify patterns visually. In addition to the standard sparsity pattern, the magnitude of the elements of the stoichiometric matrix (stoichiometric coefficients) is shown as proportional to the size of the dot.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% print a colorful spy map of the S matrix</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >spyc(S, colormap(advancedColormap(</span><span style="color: rgb(170, 4, 249);">'cobratoolbox'</span><span >)));</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% set the font size of the current figure axes</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >set(gca, </span><span style="color: rgb(170, 4, 249);">'fontsize'</span><span >, 14);</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="5CA88DA4" data-testid="output_8" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S8'><span>In the case of the </span><span style=' font-style: italic;'>E.coli core</span><span> model [2],  the biomass reaction is clearly visible (reaction number 13) due to its large amount of metabolites (dots in the column) and large coefficients (size of the dots). Also, the metabolites with large stoichiometric coefficients can be easily determined based on their dot size.</span></div><div  class = 'S1'><span style=' font-weight: bold;'>Rank</span></div><div  class = 'S1'><span>The </span><span style=' font-weight: bold;'>rank</span><span> of a stoichiometric matrix is the maximum number of linearly independent rows, and is equivalent to the number of linearly independent columns. The rank is a measurement of how many reactions and metabolites are linearly independent. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the rank of the stoichiometric matrix</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">if </span><span >ispc</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    rankS = rank(full(S))</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">else</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    rankS = getRankLUSOL(S) </span><span style="color: rgb(2, 128, 9);">% calculated using either the LUSOL solver [3]</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span style="color: rgb(14, 0, 255);">end</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>rankS = 67</div></div></div></div><div  class = 'S8'><span>The </span><span style=' font-weight: bold;'>rank deficiency</span><span> of the stoichiometric matrix is a measure of how many reactions and metabolites are not linearly dependent, and expressed as a ratio of the rank of the stoichiometric matrix to the theoretical full rank.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% calculate the rank deficiency (in percent)</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >rankDeficiencyS = (1 - rankS / min(nMets, nRxns)) * 100  </span><span style="color: rgb(2, 128, 9);">% [%]</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>rankDeficiencyS = 6.9444</div></div></div></div><div  class = 'S1'><span style=' font-weight: bold;'>Singular Values and Condition Number</span></div><div  class = 'S1'><span>A singular value decomposition of the stoichiometric matrix is the decomposition into orthonormal matrices </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">U</span><span> (of dimension </span><span style=' font-family: monospace;'>nMets</span><span> by </span><span style=' font-family: monospace;'>nMets</span><span>) and </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">V</span><span> (of dimension </span><span style=' font-family: monospace;'>nRxns</span><span> by </span><span style=' font-family: monospace;'>nRxns</span><span>), and a matrix with nonnegative diagonal elements </span><span style="font-family: STIXGeneral, STIXGeneral-webfont, serif; font-style: italic; font-weight: normal; color: rgb(0, 0, 0);">D</span><span> such that </span><span mathmlencoding="&lt;math xmlns=&quot;http://www.w3.org/1998/Math/MathML&quot; display=&quot;inline&quot;&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;normal&quot;&gt;S = UD&lt;/mi&gt;&lt;msup&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;V&lt;/mi&gt;&lt;/mrow&gt;&lt;mrow&gt;&lt;mi mathvariant=&quot;italic&quot;&gt;T&lt;/mi&gt;&lt;/mrow&gt;&lt;/msup&gt;&lt;/mrow&gt;&lt;/math&gt;" style="vertical-align:-5px"><img src="" width="68" height="19" /></span><span>.</span></div><div  class = 'S1'><span>Note that the calculation of singular values is numerically expensive, especially for large stoichiometric matrices.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% calculate the singular values</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >svVect = svds(S, rankS);</span></span></div></div></div><div  class = 'S8'><span>The </span><span style=' font-family: monospace;'>svds() </span><span>function returns the number of singular values specified in the second argument of the function. As most stoichiometric matrices are rank deficient, some singular values are zero (or within numerical tolerances). The cut-off is located at the rank of the stoichiometric matrix. </span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the vector with all singular values (including zeros)</span></span></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >svVectAll = svds(S, min(nMets, nRxns));</span></span></div></div></div><div  class = 'S8'><span>The singular values and their cut-off can be illustrated as follows:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% plot the singular values</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >figure;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% plot the singular values up to rankS</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >semilogy(linspace(1, length(svVect), length(svVect)), svVect, </span><span style="color: rgb(170, 4, 249);">'*'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% plot all singular values</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">on</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >semilogy(linspace(1, length(svVectAll), length(svVectAll)), svVectAll, </span><span style="color: rgb(170, 4, 249);">'ro'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% set the font size of the current figure axes, show a legend and minor grid axes</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >set(gca, </span><span style="color: rgb(170, 4, 249);">'fontsize'</span><span >, 14);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >legend(</span><span style="color: rgb(170, 4, 249);">'svds (up to rankS)'</span><span >, </span><span style="color: rgb(170, 4, 249);">'svds (all)'</span><span >)</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >grid </span><span style="color: rgb(170, 4, 249);">minor</span><span >;</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% set the label</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >xlabel(</span><span style="color: rgb(170, 4, 249);">'Number of the singular value'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >ylabel(</span><span style="color: rgb(170, 4, 249);">'Magnitude of the singular value'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >hold </span><span style="color: rgb(170, 4, 249);">off</span><span >;</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsFigure" uid="B6D368E7" data-testid="output_11" style="width: 450px;"><div class="figureElement"><img class="figureImage figureContainingNode" src="" style="width: 560px;"></div></div></div></div></div><div  class = 'S8'><span>The</span><span style=' font-weight: bold;'> maximum singular</span><span> value is</span><span style=' font-weight: bold;'> </span><span>the largest element on the diagonal matrix obtained from singular value decomposition. Similarly, the </span><span style=' font-weight: bold;'>minimum singular value</span><span> is the smallest element on the diagonal matrix obtained from singular value decomposition. Only singular values greater than zero (numbered from </span><span style=' font-family: monospace;'>1</span><span> to </span><span style=' font-family: monospace;'>rank(S)</span><span>) are of interest.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the maximum and minimum singular values</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >maxSingVal = svVect(1) </span><span style="color: rgb(2, 128, 9);">% first value of the vector with singular values</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>maxSingVal = 135.5764</div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span >minSingVal = svVect(rankS) </span><span style="color: rgb(2, 128, 9);">% smallest non-zero singular value</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>minSingVal = 0.1161</div></div></div></div><div  class = 'S8'><span>Alternatively, if the rank of the stoichiometric matrix </span><span style=' font-family: monospace;'>S</span><span> is not known, the built-in functions can also be used: </span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S13'><span style="white-space: pre"><span >maxSingValBuiltIn = svds(S, 1)</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>maxSingValBuiltIn = 135.5764</div></div></div><div class="inlineWrapper outputs"><div  class = 'S12'><span style="white-space: pre"><span >minSingValBuiltIn = svds(S, 1, </span><span style="color: rgb(170, 4, 249);">'smallestnz'</span><span >)</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>minSingValBuiltIn = 0.1161</div></div></div></div><div  class = 'S8'><span>The </span><span style=' font-weight: bold;'>condition number</span><span> of the stoi</span><span>chiometric matrix is a ratio of the maximum and minimum singular values. The higher this ratio, the more ill-conditioned the stoichiometric matrix is (numerical issues) and, generally, the longer the simulation time is.</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% determine the condition number</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >condNumber = maxSingVal / minSingVal</span></span></div><div  class = 'S10'><div class='variableElement' style='font-family: Menlo, Monaco, Consolas, "Courier New", monospace; font-size: 12px; '>condNumber = 1.1676e+03</div></div></div></div><div  class = 'S1'><span style=' font-weight: bold;'>Summary of model characteristics</span></div><div  class = 'S1'><span>The following numerical properties have been calculated:</span></div><ul  class = 'S14'><li  class = 'S15'><span style=' font-weight: bold;'>Number of elements</span><span>: represents the total number of entries in the stoichiometric matrix (including zero elements). This number is equivalent to the product of the number of reactions and the number of metabolites.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Number of nonzero elements</span><span>: represents the total number of nonzero entries in the stoichiometric matrix (excluding zero elements).</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Sparsity ratio</span><span>: ratio of the number of zero elements and the total number of elements. The sparser a stoichiometric matrix, the fewer metabolites participate in each reaction. The sparsity ratio is particularly useful to compare models by how many metabolites participate in each reaction.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Complementary sparsity ratio</span><span>: calculated as the difference of one and the sparsity ratio, and is the ratio of the number of nonzero elements and the total number of elements.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Average column density</span><span>: corresponds to the ratio of the number of nonzero elements in each column and the total number of metabolites. The average column density corresponds to the arithmetic average of all the column densities (sum of all the column densities divided by the number of reactions).</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Relative column density</span><span>: corresponds to the ratio of the number of nonzero elements in each column and the total number of metabolites. The relative column density corresponds to the average column density divided by the total number of metabolites (expressed in parts-per-million (ppm)).</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Rank</span><span>: the rank of a stoichiometric matrix is the maximum number of linearly independent rows and is equivalent to the number of linearly independent columns. The rank is a measurement of how many reactions and metabolites are linearly independent.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Rank deficiency</span><span>: the rank deficiency of the stoichiometric matrix is a measure of how many reactions and metabolites are linearly dependent, and expressed as the ratio of the rank of the stoichiometric matrix to the theoretical full rank.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Maximum singular value</span><span>: the largest element on the diagonal matrix obtained from singular value decomposition.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Minimum singular value</span><span>: the smallest element on the diagonal matrix obtained from singular value decomposition.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Condition number</span><span>: the condition number of the stoichiometric matrix is the ratio of the maximum and minimum singular values. The higher this ratio, the more ill-conditioned the stoichiometric matrix is (numerical issues).</span></li></ul><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span >fprintf([</span><span style="color: rgb(170, 4, 249);">' --- SUMMARY ---\n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Model file/Model name/Matrix name    %s/%s/%s\n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Size is [nMets, nRxns]               [%d, %d]\n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Number of elements:                  %d \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Number of nonzero elements:          %d \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Sparsity ratio [%%]:                  %1.2f \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Complementary sparsity ratio [%%]     %1.2f \n'</span><span >, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Average column density [ppm]:        %1.2f \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Relative column density [ppm]:       %1.2f \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Rank:                                %d \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Rank deficiency [%%]:                 %1.2f \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Maximum singular value:              %1.2f \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Minimum singular value:              %1.2f \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    </span><span style="color: rgb(170, 4, 249);">'Condition number:                    %1.2f \n'</span><span >,</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    ],</span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    modelFile, modelName, matrixName, nMets, nRxns, nElem, nNz, sparsityRatio, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >    compSparsityRatio, colDensityAv, colDensityRel, rankS, rankDeficiencyS, </span><span style="color: rgb(14, 0, 255);">...</span></span></div></div><div class="inlineWrapper outputs"><div  class = 'S9'><span style="white-space: pre"><span >    maxSingVal, minSingVal, condNumber);</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="1F6D9F7E" data-testid="output_17" data-width="420" data-height="199" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> --- SUMMARY ---
Model file/Model name/Matrix name    ecoli_core_model.mat/model/S
Size is [nMets, nRxns]               [72, 95]
Number of elements:                  6840 
Number of nonzero elements:          360 
Sparsity ratio [%]:                  94.74 
Complementary sparsity ratio [%]     5.26 
Average column density [ppm]:        3.79 
Relative column density [ppm]:       5.26 
Rank:                                67 
Rank deficiency [%]:                 6.94 
Maximum singular value:              135.58 
Minimum singular value:              0.12 
Condition number:                    1167.63 </div></div></div></div></div><div  class = 'S1'><span style=' font-weight: bold;'>Scaling</span></div><div  class = 'S1'><span>The scaling estimate is based on the order of magnitude of the ratio of the maximum and minimum scaling coefficients, which are determined such that the scaled stoichiometric matrix has entries close to unity. In order to investigate the scaling of the stoichiometric matrix and provide an estimate of the most appropriate precision of the solver to be used, the following quantities should be calculated:</span></div><ul  class = 'S14'><li  class = 'S15'><span style=' font-weight: bold;'>Estimation level: </span><span>The estimation level, defined by the parameter scltol provides a measure of how accurate the estimation is. The estimation level can be </span><span style=' font-style: italic;'>crude</span><span>, </span><span style=' font-style: italic;'>medium</span><span>, or </span><span style=' font-style: italic;'>fine</span><span>.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Size of the matrix: </span><span>The size of the matrix indicates the size of the metabolic network, and is broken down into number of metabolites and number of reactions.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Stoichiometric coefficients:</span><span> The maximum and minimum values of the stoichiometric matrix provide a range of the stoichiometric coefficients and are determined based on all elements of the stoichiometric matrix. Their ratio (and its order of magnitude) provides valuable information on the numerical difficulty to solve a linear program.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Lower bound coefficients: </span><span>The maximum and minimum values of the lower bound vector provide a range of the coefficients of the lower bound vector. Their ratio (and its order of magnitude) provides valuable information on the numerical difficulty to solve a linear program.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Upper bound coefficients: </span><span>The maximum and minimum values of the upper bound vector provide a range of the coefficients of the upper bound vector. Their ratio (and its order of magnitude) provides valuable information on the numerical difficulty to solve a linear program.</span></li><li  class = 'S15'><span style=' font-weight: bold;'>Row scaling coefficients: </span><span>The row scaling coefficients are the scaling coefficients required to scale each row closer to unity. The maximum and minimum row scaling coefficients provide a range of row scaling coefficients required to scale the stoichiometric matrix row-wise. Their ratio (and its order of magnitude) provides valuable information on the numerical difficulty to solve a linear program. </span></li><li  class = 'S15'><span style=' font-weight: bold;'>Column scaling coefficients: </span><span>The column scaling coefficients are the scaling coefficients required to scale each column closer to unity. The maximum and minimum column scaling coefficients provide a range of column scaling coefficients required to scale the stoichiometric matrix column-wise. Their ratio (and its order of magnitude) provides valuable information on the numerical difficulty to solve a linear program.</span></li></ul><div  class = 'S1'><span>The scaling properties of the stoichiometric matrix can be determined using:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S13'><span style="white-space: pre"><span >[precisionEstimate, solverRecommendation] = checkScaling(model);</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="4CDF2905" data-testid="output_18" data-width="420" data-height="619" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> ------------------------ Scaling summary report ------------------------

 Name of model:                                ecoli_core_model
 Estimation level:                             fine (scltol = 1.00)
 Name of matrix:                               S
 Size of matrix:
        * metabolites:                         72
        * reactions:                           95
 Stoichiometric coefficients:
        * Minimum (absolute non-zero value):   7.09e-02
        * Maximum (absolute non-zero value):   5.98e+01
 Lower bound coefficients:
        * Minimum (absolute non-zero value):   8.39e+00
        * Maximum (absolute non-zero value):   1.00e+03
 Upper bound coefficients:
        * Minimum (absolute non-zero value):   8.39e+00
        * Maximum (absolute non-zero value):   1.00e+03
 Row scaling coefficients:
        * Minimum:                             2.66e-01 (row #: 26)
        * Maximum:                             7.73e+00 (row #: 13)
 Column scaling coefficients:
        * Minimum:                             1.29e-01 (column #: 11)
        * Maximum:                             7.73e+00 (column #: 13)

 ---------------------------------- Ratios --------------------------------

 Ratio of stoichiometric coefficients:         8.44e+02
 Order of magnitude diff. (stoich. coeff.):    2

 Ratio of lower bounds:                        1.19e+02
 Order of magnitude diff. (lower bounds):      2

 Ratio of upper bounds:                        1.19e+02
 Order of magnitude diff. (upper bounds):      2

 Ratio of row scaling coefficients:            2.90e+01
 Order of magnitude diff. (row scaling):       1

 Ratio of column scaling coefficients:         5.98e+01
 Order of magnitude diff. (column scaling):    1

 --------------------------------------------------------------------------

 -&gt; The model is well scaled. Double precision is recommended.</div></div></div></div></div><div  class = 'S8'><span>The </span><span style=' font-family: monospace;'>precisionEstimate</span><span> yields a recommended estimate of the precision of the solver:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S13'><span style="white-space: pre"><span >precisionEstimate</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="1D478233" data-testid="output_19" data-width="420" data-height="20" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">precisionEstimate = </span>double</div></div></div></div></div></div><div  class = 'S8'><span>The solver recommendation is provided in </span><span style=' font-family: monospace;'>solverRecommendation</span><span> as a cell array that  can be used programmatically:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S13'><span style="white-space: pre"><span >solverRecommendation</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement scrollableOutput" uid="EE7079A9" data-testid="output_20" data-width="420" data-height="34" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">solverRecommendation = </span></div><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'glpk'    'gurobi'    'ibm_cplex'    'matlab'    'mosek'    'pdco'    'qpng'    'lp_solve'
</div></div></div></div></div></div><div  class = 'S8'><span>In order to see the effect of scaling, let us consider the ME model [4]:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S5'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% load the modelName structure from the modelFile</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">%as before this model is distributed and we need to make sure, that the right file is choosen.</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'><span style="white-space: pre"><span >modelFolder = getDistributedModelFolder(</span><span style="color: rgb(170, 4, 249);">'ME_matrix_GlcAer_WT.mat'</span><span >);</span></span></div></div><div class="inlineWrapper"><div  class = 'S6'></div></div><div class="inlineWrapper"><div  class = 'S7'><span style="white-space: pre"><span >modelGlcOAer_WT = readCbModel([modelFolder filesep </span><span style="color: rgb(170, 4, 249);">'ME_matrix_GlcAer_WT.mat'</span><span >], </span><span style="color: rgb(170, 4, 249);">'modelGlcOAer_WT'</span><span >);</span></span></div></div></div><div  class = 'S8'><span>The scaling can be checked as follows:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S13'><span style="white-space: pre"><span >[precisionEstimate, solverRecommendation] = checkScaling(modelGlcOAer_WT);</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsTextElement scrollableOutput" uid="482517DC" data-testid="output_21" data-width="420" data-height="633" data-hashorizontaloverflow="true" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"> ------------------------ Scaling summary report ------------------------

 Estimation level:                             fine (scltol = 1.00)
 Name of matrix:                               S
 Size of matrix:
        * metabolites:                         68298
        * reactions:                           76664
 Stoichiometric coefficients:
        * Minimum (absolute non-zero value):   5.50e-05
        * Maximum (absolute non-zero value):   8.01e+05
 Lower bound coefficients:
        * Minimum (absolute non-zero value):   1.00e+06
        * Maximum (absolute non-zero value):   1.00e+09
 Upper bound coefficients:
        * Minimum (absolute non-zero value):   5.54e+00
        * Maximum (absolute non-zero value):   1.00e+09
 Row scaling coefficients:
        * Minimum:                             1.67e-03 (row #: 24548)
        * Maximum:                             7.90e+04 (row #: 66206)
 Column scaling coefficients:
        * Minimum:                             3.47e-04 (column #: 130)
        * Maximum:                             1.13e+04 (column #: 24943)

 ---------------------------------- Ratios --------------------------------

 Ratio of stoichiometric coefficients:         1.46e+10
 Order of magnitude diff. (stoich. coeff.):    10

 Ratio of lower bounds:                        1.00e+03
 Order of magnitude diff. (lower bounds):      3

 Ratio of upper bounds:                        1.80e+08
 Order of magnitude diff. (upper bounds):      8

 Ratio of row scaling coefficients:            4.72e+07
 Order of magnitude diff. (row scaling):       7

 Ratio of column scaling coefficients:         3.27e+07
 Order of magnitude diff. (column scaling):    7

 --------------------------------------------------------------------------

 -&gt; The model has badly scaled rows and columns. Quad precision is strongly recommended.

    Set the Quad MINOS solver with: &gt;&gt; changeCobraSolver('quadMinos', 'LP')</div></div></div></div></div><div  class = 'S8'><span>In this case, the </span><span style=' font-family: monospace;'>solverRecommendation</span><span> is:</span></div><div class="CodeBlock"><div class="inlineWrapper outputs"><div  class = 'S13'><span style="white-space: pre"><span >solverRecommendation</span></span></div><div  class = 'S10'><div class="inlineElement eoOutputWrapper embeddedOutputsVariableStringElement" uid="EEAC7A55" data-testid="output_22" data-width="420" data-height="34" data-hashorizontaloverflow="false" style="width: 450px; max-height: 261px; white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div class="textElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;"><span class="variableNameElement" style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">solverRecommendation = </span></div><div style="white-space: pre; font-style: normal; color: rgb(64, 64, 64); font-size: 12px;">    'dqqMinos'    'quadMinos'
</div></div></div></div></div></div><div  class = 'S1'><span>The solver then can be set as suggested to use the quad-precision solver [5]:</span></div><div class="CodeBlock"><div class="inlineWrapper"><div  class = 'S4'><span style="white-space: pre"><span style="color: rgb(2, 128, 9);">% changeCobraSolver('dqqMinos');</span></span></div></div></div><div  class = 'S8'><span>Note that the timing for obtaining a solution using a quad-precision solver is very different than obtaining the solution using a double-precision solver. </span></div><h2  class = 'S3'><span>EXPECTED RESULTS</span></h2><div  class = 'S1'><span>The expected result is a summary table with relevant numerical characteristics and a recommendation of the precision of the solver.</span></div><h2  class = 'S3'><span>TROUBLESHOOTING</span></h2><div  class = 'S1'><span>If any numerical issues arise (e.g., infeasible solution status after performing flux balance analysis, or too long simulation times) when using a double precision solver, then a higher precision solver should be tested. For instance, a double precision solver may incorrectly report that an ill-scaled optimisation problem is infeasible although it actually might be feasible for a higher precision solver. </span></div><div  class = 'S1'><span>In some cases,  the precision recommendation might not be accurate enough and a double precision solver might lead to inaccurate results. Try a quad precision solver in order to confirm the results when in doubt.</span></div><div  class = 'S1'><span>The </span><span style=' font-family: monospace;'>checkScaling()</span><span> function may be used on all operating systems, but the </span><span style=' font-family: monospace;'>'dqqMinos'</span><span> interface is only available on UNIX operating systems. If the </span><span style=' font-family: monospace;'>'dqqMinos'</span><span> interface is not working as intended, the binaries might not be compatible (raise an issue if they are not by providing the output of </span><span style=' font-family: monospace;'>generateSystemConfigReport</span><span>). Make sure that all relevant system requirements are satisfied before trying to use the </span><span style=' font-family: monospace;'>'dqqMinos'</span><span> solver interface.</span></div><div  class = 'S1'><span>In case the </span><span style=' font-family: monospace;'>'dqqMinos' </span><span>interface reports an error when trying to solve the linear program, there might be an issue with the model itself.</span></div><h2  class = 'S3'><span>References</span></h2><div  class = 'S1'><span>[1] </span><a href = "http://www.nature.com/nbt/journal/v31/n5/full/nbt.2488.html"><span style=' text-decoration: underline;'>Thiele et al., A community-driven global reconstruction of human metabolism, Nat Biotech, 2013.</span></a></div><div  class = 'S1'><span>[2] Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia coli Metabolic Model as an Educational Guide by Orth, Fleming, and Palsson (2010)</span></div><div  class = 'S1'><span>[3] P. E. Gill, W. Murray, M. A. Saunders and M. H. Wright (1987). Maintaining LU factors of a general sparse matrix, Linear Algebra and its Applications 88/89, 239-270.</span></div><div  class = 'S1'><span>[4] Multiscale modeling of metabolism and macromolecular synthesis in E. coli and its application to the evolution of codon usage, Thiele et al., PLoS One, 7(9):e45635 (2012).</span></div><div  class = 'S1'><span>[5] D. Ma, L. Yang, R. M. T. Fleming, I. Thiele, B. O. Palsson and M. A. Saunders, Reliable and efficient solution of genome-scale models of Metabolism and macromolecular Expression, Scientific Reports 7, 40863; doi: \url{10.1038/srep40863} (2017). </span><a href = "http://rdcu.be/oCpn"><span style=' text-decoration: underline;'>http://rdcu.be/oCpn</span></a><span>.</span></div>
<br>
<!-- 
##### SOURCE BEGIN #####
%% *Numerical properties of a reconstruction*
% *Authors: Laurent Heirendt, Ronan M.T. Fleming, Luxembourg Centre for Systems 
% Biomedicine*
% 
% *Reviewers: Sylvain Arreckx, Thomas Pfau, and Catherine Fleming,  Luxembourg 
% Centre for Systems Biomedicine*
%% Introduction
% During this tutorial, you will learn how to determine and explore the numerical 
% properties of a stoichiometric matrix. The numerical properties are key to analyzing 
% the metabolic reconstruction at hand, to select the appropriate solver, or to 
% determine incoherencies in the network. 
% 
% You will also be able to learn more about the definitions of the various numerical 
% characteristics. This tutorial is particularly useful when you have a multi-scale 
% model and are facing numerical issues when performing flux balance analysis 
% or any other variants of FBA.
%% EQUIPMENT SETUP
%% *Initialize the COBRA Toolbox.*
% Please ensure that The COBRA Toolbox has been properly installed, and initialized 
% using the |initCobraToolbox| function.

initCobraToolbox(false) % false, as we don't want to update
%% PROCEDURE 
% TIMING: 5 seconds - several hours (depending on the model size)
% 
% *Define the name of the model*
% 
% Throughout this tutorial, we will use the _E.coli core_ model [2]. It is generally 
% good practice to define the name of the file that contains the model, the name 
% of the model structure, and the name of the stoichiometric matrix, as separate 
% variables. Therefore, we propose that within the |modelFile|, there is a structure 
% named |modelName| with a field |matrixName| that contains the stoichiometric 
% matrix |S| (or |A|).

% If this is a distributed model, get the folder for the model. Commonly if you use your own model, this is unnecessary
% As the file would be in the current folder. But for this tutorial we need to make sure, that the right model is used, 
% and that no other model with the same name is above it in the path.
modelFolder = getDistributedModelFolder('ecoli_core_model.mat');

% define the filename of the model
modelFile = [modelFolder filesep 'ecoli_core_model.mat'];

% define the name of model structure
modelName = 'model';

% define the fieldname of the stoichiometric matrix
matrixName = 'S';
%% 
% *Load the stoichiometric matrix*
% 
% In order to use the model, we need to read the |modelFile| that contains a 
% COBRA model structure  |modelName|:

% load the modelName structure from the modelFile
model = readCbModel(modelFile, 'modelName','model');
%% 
% Some models contain stoichiometric matrices with a different name (commonly 
% coupled models). By default, the stoichiometric matrix is denoted |S|.

% select the matrix
S = model.S;
if isfield(model, matrixName) == 1 && strcmp(matrixName, 'A') == 1
    S = model.A;
end
%% 
% *Basic numerical characteristics*
% 
% The *number of elements* represents the total number of entries in the stoichiometric 
% matrix (including zero elements). This number is equivalent to the product of 
% the number of reactions and the number of metabolites.
% 
% The number of rows represents the *number of metabolites* in the metabolic 
% network. The number of columns corresponds to the *number of biochemical reactions* 
% in the network.

% determine the number of reactions and metabolites in S
[nMets, nRxns] = size(S)
% determine the number of elements in S
nElem = numel(S)  % Nmets * Nrxns
%% 
% The total number of nonzero elements corresponds to the total number of nonzero 
% entries in the stoichiometric matrix (excluding zero elements).

% determine the number of nonzero elements in S
nNz = nnz(S)
%% 
% *Sparsity and Density*
% 
% The *sparsity ratio* corresponds to a ratio of the number of zero elements 
% and the total number of elements. The sparser a stoichiometric matrix, the fewer 
% metabolites participate in each reaction. The sparsity ratio is particularly 
% useful to compare models by how many metabolites participate in each reaction.

% determine the sparsity ratio of S (in percent)
sparsityRatio = (1 - nNz / nElem) * 100.0  % [%]
%% 
% Similarly, the *complementary sparsity ratio* is calculated as the difference 
% of 100 and the sparsity ratio expressed in percent, and therefore, is a ratio 
% of the number of nonzero elements and the total number of elements.

% determine the complementary sparsity ratio (in percent)
compSparsityRatio = 100.0 - sparsityRatio  % [%]
%% 
% The *average column density* corresponds to a ratio of the number of nonzero 
% elements in each column (i.e. reaction) and the total number of metabolites. 
% The average column density corresponds to the arithmetic average of all the 
% column densities (sum of all the reaction densities divided by the number of 
% reactions).

% add the number of non-zeros in each column (reaction)
colDensityAv = 0;
for j = 1:nRxns
    colDensityAv = colDensityAv + nnz(S(:, j));
end

% calculate the arithmetic average number of non-zeros in each column
colDensityAv = colDensityAv / nRxns   % [-]
%% 
% The average column density provides a measure of how many stoichiometric coefficients 
% participate in each biochemical reaction on average.
% 
% The *relative column density* corresponds to the ratio of the number of nonzero 
% elements in each column and the total number of metabolites. The relative column 
% density corresponds to the average column density divided by the total number 
% of metabolites (expressed in percent). The relative column density may also 
% be expressed as parts-per-million [ppm] for large-scale or huge-scale models. 

% determine the density proportional to the length of the column
colDensityRel = colDensityAv / nMets * 100  % [%]
%% 
% The relative column density indicates how many metabolites are being used 
% on average in each reaction relative to the total number of metabolites in the 
% metabolic network.
%% 
% *Sparsity Pattern (spy plot)*
% 
% The visualisation of the sparsity pattern is useful to explore the matrix, 
% spot inconsistencies, or identify patterns visually. In addition to the standard 
% sparsity pattern, the magnitude of the elements of the stoichiometric matrix 
% (stoichiometric coefficients) is shown as proportional to the size of the dot.

% print a colorful spy map of the S matrix
spyc(S, colormap(advancedColormap('cobratoolbox')));

% set the font size of the current figure axes
set(gca, 'fontsize', 14);
%% 
% In the case of the _E.coli core_ model [2],  the biomass reaction is clearly 
% visible (reaction number 13) due to its large amount of metabolites (dots in 
% the column) and large coefficients (size of the dots). Also, the metabolites 
% with large stoichiometric coefficients can be easily determined based on their 
% dot size.
%% 
% *Rank*
% 
% The *rank* of a stoichiometric matrix is the maximum number of linearly independent 
% rows, and is equivalent to the number of linearly independent columns. The rank 
% is a measurement of how many reactions and metabolites are linearly independent. 

% determine the rank of the stoichiometric matrix
if ispc
    rankS = rank(full(S))
else
    rankS = getRankLUSOL(S) % calculated using either the LUSOL solver [3]
end
%% 
% The *rank deficiency* of the stoichiometric matrix is a measure of how many 
% reactions and metabolites are not linearly dependent, and expressed as a ratio 
% of the rank of the stoichiometric matrix to the theoretical full rank.

% calculate the rank deficiency (in percent)
rankDeficiencyS = (1 - rankS / min(nMets, nRxns)) * 100  % [%]
%% 
% *Singular Values and Condition Number*
% 
% A singular value decomposition of the stoichiometric matrix is the decomposition 
% into orthonormal matrices $U$ (of dimension |nMets| by |nMets|) and $V$ (of 
% dimension |nRxns| by |nRxns|), and a matrix with nonnegative diagonal elements 
% $D$ such that $\textrm{S}\;\textrm{=}\;\textrm{UD}V^T$.
% 
% Note that the calculation of singular values is numerically expensive, especially 
% for large stoichiometric matrices.

% calculate the singular values
svVect = svds(S, rankS);
%% 
% The |svds()| function returns the number of singular values specified in the 
% second argument of the function. As most stoichiometric matrices are rank deficient, 
% some singular values are zero (or within numerical tolerances). The cut-off 
% is located at the rank of the stoichiometric matrix. 

% determine the vector with all singular values (including zeros)
svVectAll = svds(S, min(nMets, nRxns));
%% 
% The singular values and their cut-off can be illustrated as follows:

% plot the singular values
figure;

% plot the singular values up to rankS
semilogy(linspace(1, length(svVect), length(svVect)), svVect, '*');

% plot all singular values
hold on;
semilogy(linspace(1, length(svVectAll), length(svVectAll)), svVectAll, 'ro');

% set the font size of the current figure axes, show a legend and minor grid axes
set(gca, 'fontsize', 14);
legend('svds (up to rankS)', 'svds (all)')
grid minor;

% set the label
xlabel('Number of the singular value');
ylabel('Magnitude of the singular value');

hold off;
%% 
% The *maximum singular* value is the largest element on the diagonal matrix 
% obtained from singular value decomposition. Similarly, the *minimum singular 
% value* is the smallest element on the diagonal matrix obtained from singular 
% value decomposition. Only singular values greater than zero (numbered from |1| 
% to |rank(S)|) are of interest.

% determine the maximum and minimum singular values
maxSingVal = svVect(1) % first value of the vector with singular values
minSingVal = svVect(rankS) % smallest non-zero singular value
%% 
% Alternatively, if the rank of the stoichiometric matrix |S| is not known, 
% the built-in functions can also be used: 

maxSingValBuiltIn = svds(S, 1)
minSingValBuiltIn = svds(S, 1, 'smallestnz')
%% 
% The *condition number* of the stoichiometric matrix is a ratio of the maximum 
% and minimum singular values. The higher this ratio, the more ill-conditioned 
% the stoichiometric matrix is (numerical issues) and, generally, the longer the 
% simulation time is.

% determine the condition number
condNumber = maxSingVal / minSingVal
%% 
% *Summary of model characteristics*
% 
% The following numerical properties have been calculated:
%% 
% * *Number of elements*: represents the total number of entries in the stoichiometric 
% matrix (including zero elements). This number is equivalent to the product of 
% the number of reactions and the number of metabolites.
% * *Number of nonzero elements*: represents the total number of nonzero entries 
% in the stoichiometric matrix (excluding zero elements).
% * *Sparsity ratio*: ratio of the number of zero elements and the total number 
% of elements. The sparser a stoichiometric matrix, the fewer metabolites participate 
% in each reaction. The sparsity ratio is particularly useful to compare models 
% by how many metabolites participate in each reaction.
% * *Complementary sparsity ratio*: calculated as the difference of one and 
% the sparsity ratio, and is the ratio of the number of nonzero elements and the 
% total number of elements.
% * *Average column density*: corresponds to the ratio of the number of nonzero 
% elements in each column and the total number of metabolites. The average column 
% density corresponds to the arithmetic average of all the column densities (sum 
% of all the column densities divided by the number of reactions).
% * *Relative column density*: corresponds to the ratio of the number of nonzero 
% elements in each column and the total number of metabolites. The relative column 
% density corresponds to the average column density divided by the total number 
% of metabolites (expressed in parts-per-million (ppm)).
% * *Rank*: the rank of a stoichiometric matrix is the maximum number of linearly 
% independent rows and is equivalent to the number of linearly independent columns. 
% The rank is a measurement of how many reactions and metabolites are linearly 
% independent.
% * *Rank deficiency*: the rank deficiency of the stoichiometric matrix is a 
% measure of how many reactions and metabolites are linearly dependent, and expressed 
% as the ratio of the rank of the stoichiometric matrix to the theoretical full 
% rank.
% * *Maximum singular value*: the largest element on the diagonal matrix obtained 
% from singular value decomposition.
% * *Minimum singular value*: the smallest element on the diagonal matrix obtained 
% from singular value decomposition.
% * *Condition number*: the condition number of the stoichiometric matrix is 
% the ratio of the maximum and minimum singular values. The higher this ratio, 
% the more ill-conditioned the stoichiometric matrix is (numerical issues).

fprintf([' REPLACE_WITH_DASH_DASH- SUMMARY REPLACE_WITH_DASH_DASH-\n',...
    'Model file/Model name/Matrix name    %s/%s/%s\n',...
    'Size is [nMets, nRxns]               [%d, %d]\n',...
    'Number of elements:                  %d \n',...
    'Number of nonzero elements:          %d \n',...
    'Sparsity ratio [%%]:                  %1.2f \n',...
    'Complementary sparsity ratio [%%]     %1.2f \n', ...
    'Average column density [ppm]:        %1.2f \n',...
    'Relative column density [ppm]:       %1.2f \n',...
    'Rank:                                %d \n',...
    'Rank deficiency [%%]:                 %1.2f \n',...
    'Maximum singular value:              %1.2f \n',...
    'Minimum singular value:              %1.2f \n',...
    'Condition number:                    %1.2f \n',...
    ],...
    modelFile, modelName, matrixName, nMets, nRxns, nElem, nNz, sparsityRatio, ...
    compSparsityRatio, colDensityAv, colDensityRel, rankS, rankDeficiencyS, ...
    maxSingVal, minSingVal, condNumber);
%% 
% *Scaling*
% 
% The scaling estimate is based on the order of magnitude of the ratio of the 
% maximum and minimum scaling coefficients, which are determined such that the 
% scaled stoichiometric matrix has entries close to unity. In order to investigate 
% the scaling of the stoichiometric matrix and provide an estimate of the most 
% appropriate precision of the solver to be used, the following quantities should 
% be calculated:
%% 
% * *Estimation level:* The estimation level, defined by the parameter scltol 
% provides a measure of how accurate the estimation is. The estimation level can 
% be _crude_, _medium_, or _fine_.
% * *Size of the matrix:* The size of the matrix indicates the size of the metabolic 
% network, and is broken down into number of metabolites and number of reactions.
% * *Stoichiometric coefficients:* The maximum and minimum values of the stoichiometric 
% matrix provide a range of the stoichiometric coefficients and are determined 
% based on all elements of the stoichiometric matrix. Their ratio (and its order 
% of magnitude) provides valuable information on the numerical difficulty to solve 
% a linear program.
% * *Lower bound coefficients:* The maximum and minimum values of the lower 
% bound vector provide a range of the coefficients of the lower bound vector. 
% Their ratio (and its order of magnitude) provides valuable information on the 
% numerical difficulty to solve a linear program.
% * *Upper bound coefficients:* The maximum and minimum values of the upper 
% bound vector provide a range of the coefficients of the upper bound vector. 
% Their ratio (and its order of magnitude) provides valuable information on the 
% numerical difficulty to solve a linear program.
% * *Row scaling coefficients:* The row scaling coefficients are the scaling 
% coefficients required to scale each row closer to unity. The maximum and minimum 
% row scaling coefficients provide a range of row scaling coefficients required 
% to scale the stoichiometric matrix row-wise. Their ratio (and its order of magnitude) 
% provides valuable information on the numerical difficulty to solve a linear 
% program. 
% * *Column scaling coefficients:* The column scaling coefficients are the scaling 
% coefficients required to scale each column closer to unity. The maximum and 
% minimum column scaling coefficients provide a range of column scaling coefficients 
% required to scale the stoichiometric matrix column-wise. Their ratio (and its 
% order of magnitude) provides valuable information on the numerical difficulty 
% to solve a linear program.
%% 
% The scaling properties of the stoichiometric matrix can be determined using:

[precisionEstimate, solverRecommendation] = checkScaling(model);
%% 
% The |precisionEstimate| yields a recommended estimate of the precision of 
% the solver:

precisionEstimate
%% 
% The solver recommendation is provided in |solverRecommendation| as a cell 
% array that  can be used programmatically:

solverRecommendation
%% 
% In order to see the effect of scaling, let us consider the ME model [4]:

% load the modelName structure from the modelFile
%as before this model is distributed and we need to make sure, that the right file is choosen.
modelFolder = getDistributedModelFolder('ME_matrix_GlcAer_WT.mat');

modelGlcOAer_WT = readCbModel([modelFolder filesep 'ME_matrix_GlcAer_WT.mat'], 'modelGlcOAer_WT');
%% 
% The scaling can be checked as follows:

[precisionEstimate, solverRecommendation] = checkScaling(modelGlcOAer_WT);
%% 
% In this case, the |solverRecommendation| is:

solverRecommendation
%% 
% The solver then can be set as suggested to use the quad-precision solver [5]:

% changeCobraSolver('dqqMinos');
%% 
% Note that the timing for obtaining a solution using a quad-precision solver 
% is very different than obtaining the solution using a double-precision solver. 
%% EXPECTED RESULTS
% The expected result is a summary table with relevant numerical characteristics 
% and a recommendation of the precision of the solver.
%% TROUBLESHOOTING
% If any numerical issues arise (e.g., infeasible solution status after performing 
% flux balance analysis, or too long simulation times) when using a double precision 
% solver, then a higher precision solver should be tested. For instance, a double 
% precision solver may incorrectly report that an ill-scaled optimisation problem 
% is infeasible although it actually might be feasible for a higher precision 
% solver. 
% 
% In some cases,  the precision recommendation might not be accurate enough 
% and a double precision solver might lead to inaccurate results. Try a quad precision 
% solver in order to confirm the results when in doubt.
% 
% The |checkScaling()| function may be used on all operating systems, but the 
% |'dqqMinos'| interface is only available on UNIX operating systems. If the |'dqqMinos'| 
% interface is not working as intended, the binaries might not be compatible (raise 
% an issue if they are not by providing the output of |generateSystemConfigReport|). 
% Make sure that all relevant system requirements are satisfied before trying 
% to use the |'dqqMinos'| solver interface.
% 
% In case the |'dqqMinos'| interface reports an error when trying to solve the 
% linear program, there might be an issue with the model itself.
%% References
% [1] <http://www.nature.com/nbt/journal/v31/n5/full/nbt.2488.html Thiele et 
% al., A community-driven global reconstruction of human metabolism, Nat Biotech, 
% 2013.>
% 
% [2] Reconstruction and Use of Microbial Metabolic Networks: the Core Escherichia 
% coli Metabolic Model as an Educational Guide by Orth, Fleming, and Palsson (2010)
% 
% [3] P. E. Gill, W. Murray, M. A. Saunders and M. H. Wright (1987). Maintaining 
% LU factors of a general sparse matrix, Linear Algebra and its Applications 88/89, 
% 239-270.
% 
% [4] Multiscale modeling of metabolism and macromolecular synthesis in E. coli 
% and its application to the evolution of codon usage, Thiele et al., PLoS One, 
% 7(9):e45635 (2012).
% 
% [5] D. Ma, L. Yang, R. M. T. Fleming, I. Thiele, B. O. Palsson and M. A. Saunders, 
% Reliable and efficient solution of genome-scale models of Metabolism and macromolecular 
% Expression, Scientific Reports 7, 40863; doi: \url{10.1038/srep40863} (2017). 
% <http://rdcu.be/oCpn http://rdcu.be/oCpn>.
##### SOURCE END #####
-->
</div></body></html>